%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define the ordering in the VAR
inflpos = 1;
y1pos = 2;
ysprpos = 3;
gdppos = 4;
divgrpos = 5;

pdpos = 6;
dtpos = 7;
coint_tax = 8;
dgpos = 9;
coint_spending = 10;

% VAR elements
infl = inflation - pi0;
x = rpcgdpgr - x0;
yield1 = ynom1 - y0nom_1;
nomyspr = yieldspr - yspr0;
dg = deltalogg - g0;
dt = deltalogtau - tau0;
pd = pdm - A0m;
% divgr = divgrm - mu_m;

tts = tts';
gts = gts';

% demean all other variables
X2(:, inflpos) = [infl(date <= 1860) - mean(infl(date <= 1860)); infl(date > 1860) - mean(infl(date > 1860))];
X2(:, y1pos) = [yield1(date <= 1860) - mean(yield1(date <= 1860)); yield1(date > 1860) - mean(yield1(date > 1860))];
X2(:, ysprpos) = [nomyspr(date <= 1860) - mean(nomyspr(date <= 1860)); nomyspr(date > 1860) - mean(nomyspr(date > 1860))];
X2(:, gdppos) = [x(date <= 1860) - mean(x(date <= 1860)); x(date > 1860) - mean(x(date > 1860))];
X2(:, divgrpos) = [divgr(date <= 1860) - mean(divgr(date <= 1860)); divgr(date > 1860) - mean(divgr(date > 1860))];
X2(:, pdpos) = [pd(date <= 1860) - mean(pd(date <= 1860)); pd(date > 1860) - mean(pd(date > 1860))];
X2(:, dtpos) = [dt(date <= 1860) - mean(dt(date <= 1860)); dt(date > 1860) - mean(dt(date > 1860))];
X2(:, coint_tax) = [log(tts(date <= 1860)) - mean(log(tts(date <= 1860))); log(tts(date > 1860)) - mean(log(tts(date > 1860)))];
X2(:, dgpos) = [dg(date <= 1860) - mean(dg(date <= 1860)); dg(date > 1860) - mean(dg(date > 1860))];
X2(:, coint_spending) = [log(gts(date <= 1860)) - mean(log(gts(date <= 1860))); log(gts(date > 1860)) - mean(log(gts(date > 1860)))]';

X2(:, y1pos) = yield1;
X2(:, ysprpos) = nomyspr;
X2(:, gdppos) = x;
X2(:, divgrpos) = divgrm - mean(divgrm);
X2(:, pdpos) = pd;
X2(:, dtpos) = dt;
X2(:, coint_tax) = log(tts) - mean(log(tts));
X2(:, dgpos) = dg;
X2(:, coint_spending) = log(gts) - mean(log(gts));

N = cols(X2);
I = eye(N);
I_pi = I(:, inflpos);
I_gdp = I(:, gdppos);
I_y1 = I(:, y1pos);
I_yspr = I(:, ysprpos);
I_pdm = I(:, pdpos);
I_divgrm = I(:, divgrpos);
I_dt = I(:, dtpos);
I_dg = I(:, dgpos);
I_coint_tax = I(:, coint_tax);
I_coint_spending = I(:, coint_spending);

z_var = X2(2:end, :);
z_varlag = X2(1:end - 1, :);
Y = z_var;
X = z_varlag;

% actual, non-detrended series of tax and spending
X2t = X2;
X2t(:, coint_tax) = log(taxrevgdp) - mean(log(taxrevgdp));
X2t(:, coint_spending) = log(spendgdp) - mean(log(spendgdp));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate Psi and Sig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psi = zeros(N, N);
Psi_se = zeros(N, N);

for i = [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos dgpos]
    regr = ols(Y(:, i), [ones(T - 1, 1), X(:, :)]);
    c(i) = regr.beta(1);
    c_se(i) = regr.bstd(1);
    Psi(i, :) = regr.beta(2:end)';
    Psi_se(i, :) = regr.bstd(2:end)';
    R2(i) = regr.rsqr;
    R2bar(i) = regr.rbar;
    eps(:, i) = Y(:, i) - c(i) - X * Psi(i, :)';

    clear regr;
end

Psi(coint_tax, :) = Psi(dtpos, :);
Psi(coint_spending, :) = Psi(dgpos, :);

Psi(coint_tax, coint_tax) = Psi(coint_tax, coint_tax) + 1;
Psi(coint_spending, coint_spending) = Psi(coint_spending, coint_spending) + 1;

Psi_se(coint_tax, :) = Psi_se(dtpos, :);
Psi_se(coint_spending, :) = Psi_se(dgpos, :);

Sigma = cov(eps(:, [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos dgpos]));
tmp = chol(Sigma, 'lower');
Sig = zeros(N, N);
Sig([inflpos y1pos ysprpos gdppos divgrpos], [inflpos y1pos ysprpos gdppos divgrpos]) = ...
    tmp(1:5, 1:5);
Sig(pdpos, [inflpos y1pos ysprpos gdppos divgrpos pdpos]) = tmp(6, 1:6);
Sig(dtpos, [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos]) = tmp(7, 1:7);
Sig(dgpos, [inflpos y1pos ysprpos gdppos divgrpos pdpos dtpos dgpos]) = tmp(8, 1:8);

Sig(coint_tax, :) = Sig(dtpos, :);
Sig(coint_spending, :) = Sig(dgpos, :);

eps = [eps(:, [inflpos y1pos ysprpos gdppos divgrpos]), ...
                   eps(:, pdpos), eps(:, dtpos), eps(:, dtpos), ...
                   eps(:, dgpos), eps(:, dgpos)];
Sigma = Sig * Sig';

tstat = Psi ./ Psi_se;

max(abs(eig(Psi)))

eps2 = eps;

tau0 = 0;
g0 = 0;
